#include "FEMSpace.h"
#include <typeinfo>
#include <math.h>
#include <ctime>
#include "Error.h"
#define pi 4.0*atan(1.0)
double f(double *p)
{
    return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
}

double bc(double *p)
{
    return sin(pi * p[0]) * sin(pi * p[1]);
}

double f2(double *p)
{
    return 2 * pi*  pi* cos(pi * p[0]) * cos(pi * p[1]);
}

double bc2(double *p)
{
    return cos(pi * p[0]) * cos(pi * p[1]);
}
double a(double *p)
{
    return 1;
}
double g1(double *p)
{
    return exp(p[0] - 1.0);
}
double g2(double *p)
{
    return exp(1.0 + p[1]);
}
double g3(double *p)
{
    return exp(p[0] + 1.0);
}
double g4(double *p)
{
    return exp(-1.0 + p[1]);
}
double u(double *p)
{
    return exp(p[0] + p[1]);
}
double fe(double *p)
{
    return -2 * exp(p[0] + p[1]);
}
int main()
{
    //g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    RectangleDomain* r = new RectangleDomain({{-1,-1},{1,-1},{1,1},{-1,1}});
    std::vector<Boundary<2> > B = r->boundary();
    int segmentx = 16;
    int segmenty = segmentx;
    Mesh<2>* m = new Q1Mesh(r,{segmentx,segmenty});
    Element<2>* e = new Quadrilateral_1_Element();
    Equation<2>* equ = new PossionEquation<2>();
    equ ->SetRightHandsTermFunction(fe);
    BoundaryFunction<2> * bf1 = new Dirichlet<2>(g1,B[0]);
    BoundaryFunction<2> * bf2 = new Dirichlet<2>(g2,B[1]);
    BoundaryFunction<2> * bf3 = new Dirichlet<2>(g3,B[2]);
    BoundaryFunction<2> * bf4 = new Dirichlet<2>(g4,B[3]);
    BoundaryCondition<2> bfc;
    bfc.add(bf1);
    bfc.add(bf2);
    bfc.add(bf3);
    bfc.add(bf4);
    Possion_2D possionproblem(m,e,equ,bfc);
    possionproblem.AssembleStiffMatrix();
    possionproblem.AssembleRightHandsTerm();
    possionproblem.DealWithBoubdaryCondition();
    possionproblem.Solve();
    Eigen::VectorXd residual = possionproblem.A() * possionproblem.solution() - possionproblem.Rhs();
    double residual_l2 = residual.squaredNorm();
    double residual_max = residual.lpNorm<Eigen::Infinity>();
    std::cout<<"l2 residual is :" << residual_l2 << " ,residual max norm is :" << residual_max<<std::endl;
    Eigen::VectorXd Truesolution = Eigen::MatrixXd::Zero(m->n_dofs(),1);
    for(int j = 0; j<=segmenty;j++)
    {
        for(int i = 0;i <= segmentx;i++)
        {
            double point[2] = {-1 + i * m->x_h(),-1 + j * m->y_h()};
            Truesolution[i + (segmentx+1) * j] = u(point);
        }
    }
    std::cout << "done" << std::endl;
    Eigen::VectorXd error = Truesolution - possionproblem.solution();
    //std::cout << "error is :" << error <<std::endl;
    double error_l2 = error.norm() / segmentx;
    double error_max = error.lpNorm<Eigen::Infinity>();
    std::cout<<" Error l2 norm is :" << error_l2 << " ,error max norm is :" << error_max<<std::endl;
}